The Cancer Genome Atlas Program (TCGA) è un programma americano che ha caratterizzato 33 tipi di tumori e messo a disposizione della comunità scientifica i dati risultanti. L’obiettivo principale del progetto era creare un catalogo completo delle mutazioni genetiche e delle alterazioni molecolari responsabili dei vari tipi di cancro. Tra questi è incluso BRCA.
Linkedomics è un portale web che ha integrato i dati del TCGA e di un altro programma, e mette a disposizione dati già pre-processati.
Per svolgere il seguente tutorial è necessario installare:
il software R disponibile per il proprio sistema operativo CRAN IT Mirror, e
l’IDE RStudio distribuito da posit.
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("TCGAbiolinks")
library("TCGAbiolinks")
require(reader)
require(dplyr)
require(janitor)
require(ggplot2)
require(survival)
require(survminer)
require(glmnet)
Il pacchetto R/Bioconductor TCGA Biolinks permette di
scaricare direttamente in R i dataset presenti in Linkedomics.
# Query per i dati clinici
TCGA_BRCA_clinical = getLinkedOmicsData(
project = "TCGA-BRCA",
dataset = "Clinical"
)
dim(TCGA_BRCA_clinical)
## [1] 20 1098
# Visualizziamo i nomi delle colonne (attributi)
TCGA_BRCA_clinical$attrib_name
## [1] "years_to_birth" "Tumor_purity"
## [3] "pathologic_stage" "pathology_T_stage"
## [5] "pathology_N_stage" "pathology_M_stage"
## [7] "histological_type" "number_of_lymph_nodes"
## [9] "PAM50" "ER.Status"
## [11] "PR.Status" "HER2.Status"
## [13] "gender" "radiation_therapy"
## [15] "race" "ethnicity"
## [17] "Median_overall_survival" "overall_survival"
## [19] "status" "overallsurvival"
# Visualizziamo i dati di sopravvivenza (righe 18-20, per semplificare colonne 1-5)
TCGA_BRCA_clinical[18:20, 1:5]
## # A tibble: 3 × 5
## attrib_name TCGA.5L.AAT0 TCGA.5L.AAT1 TCGA.A1.A0SP TCGA.A2.A04V
## <chr> <chr> <chr> <chr> <chr>
## 1 overall_survival 1477 1471 584 1920
## 2 status 0 0 0 1
## 3 overallsurvival 1477,0 1471,0 584,0 1920,1
Nota. In caso di “overall_survival”, l’evento considerato nell’ analisi di sopravvivenza è la morte del paziente.
# Trasponiamo e puliamo i dati per avere un dataframe pazienti x dati clinici
TCGA_BRCA_clinical = t(TCGA_BRCA_clinical) %>%
as.data.frame() %>%
janitor::row_to_names(row_number = 1) %>%
mutate(across(everything(), ~ type.convert(as.character(.), as.is = TRUE)))
dim(TCGA_BRCA_clinical)
## [1] 1097 20
# Vediamo quanti NA abbiamo nei dati di sopravvvenza
print(paste("NA tempi sopravvivenza:", sum(is.na(TCGA_BRCA_clinical$overall_survival))))
## [1] "NA tempi sopravvivenza: 42"
print(paste("NA eventi:", sum(is.na(TCGA_BRCA_clinical$status))))
## [1] "NA eventi: 42"
# Rimuoviamo i pazienti con dati di sopravvivenza mancanti
TCGA_BRCA_clinical = TCGA_BRCA_clinical[
!is.na(TCGA_BRCA_clinical$overall_survival) &
!is.na(TCGA_BRCA_clinical$status),]
dim(TCGA_BRCA_clinical)
## [1] 1055 20
Il censoring plot ci mostra un summary dei dati di sopravvivenza, possiamo visualizzare i tempi (ascisse) e la presenza dell’evento o del censoring.
Usiamo solo i primi 30 pazienti per un plot più chiaro.
clin = TCGA_BRCA_clinical[1:30,]
clin = clin %>% arrange(overall_survival) # ordiniamo per tempi
ID = 1:nrow(clin)
clin %>%
ggplot(
aes(y = ID, x = overall_survival, shape = factor(status), color = factor(status))
) +
geom_segment(aes(x = overall_survival, y = ID, xend = 0, yend = ID), linewidth = 2) +
geom_point(size=3) +
labs(x = "Days", y = "Patients") +
scale_shape_discrete(name = "Status", labels = c("Censored","Dead")) +
scale_color_manual(
name = "Status",
values = c("0" = "skyblue1", "1" = "red"),
labels = c("0" = "Censored", "1" = "Dead")
) +
scale_x_continuous(
breaks = seq(0, max(clin$overall_survival, na.rm = TRUE), by = 180),
limits = c(0, max(clin$overall_survival, na.rm = TRUE))
) +
theme(
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
)
La funzione Surv() dal pacchetto survival
combina i due vettori:
# Creare l'oggetto Surv
surv_object = Surv(
time = TCGA_BRCA_clinical$overall_survival,
event = TCGA_BRCA_clinical$status
)
# Visualizziamo i primi 20
head(surv_object, 20)
## [1] 1477+ 1471+ 584+ 1920 1099+ 2695+ 595+ 3276+ 1165+ 718+ 954+ 738+
## [13] 724+ 661+ 679+ 343+ 364+ 336+ 1614+ 883
L’output (es. 1477+) indica un paziente seguito per 1477 giorni e censurato (vivo all’ultimo follow-up). 883 indica un paziente deceduto (evento) al giorno 883.
La curva di Kaplan-Meier è un grafico che mostra la probabilità di sopravvivenza stimata sull’asse Y in funzione del tempo sull’asse X.
Funzione a gradini: La curva è piatta tra un evento e l’altro e scende verticalmente solo nei momenti in cui si verificano uno o più eventi. L’altezza del “gradino” verso il basso dipende dal numero di eventi rispetto al numero di persone a rischio in quel momento.
Inizio al 100%: Al tempo 0, la probabilità di sopravvivenza è 1 (o 100%), poiché nessun evento si è ancora verificato.
Indicatori di censura: I tempi in cui i soggetti vengono censurati sono indicati sulla curva con piccoli segni per mostrare dove l’informazione è stata persa senza che si verificasse un evento.
fit_km_global = survfit(
surv_object ~ 1,
data = TCGA_BRCA_clinical
)
ggsurvplot(
fit_km_global,
data = TCGA_BRCA_clinical,
conf.int = TRUE,
risk.table = TRUE,
title = "Curva di Sopravvivenza Globale (TCGA-BRCA)",
xlab = "Tempo (Giorni)"
)
Creare curve separate per diversi gruppi ci permette di confrontare
l’andamento della sopravvivenza rispetto una condizione (trattamento,
stadiazione, sano vs malato).
Oltre ad una valutazione visiva, il log-rank test ci
permette di stabilire se la differenza tra le curve è significativa.
La curva che rimane costantemente sopra l’ altra rappresenta il gruppo con la sopravvivenza migliore (cioè una minore probabilità di sperimentare l’evento in ogni dato momento).
Stratifichiamo i pazienti rispetto il dato clinico radioterapia (yes/no).
#Rimuoviamo i pazienti con il dato clinico mancante
print(paste("NA radiation_therapy:", sum(is.na(TCGA_BRCA_clinical$radiation_therapy))))
## [1] "NA radiation_therapy: 79"
TCGA_BRCA_clinical_clean = TCGA_BRCA_clinical[!is.na(TCGA_BRCA_clinical$radiation_therapy), ]
fit_km_radiation_therapy = survfit(
Surv(overall_survival, status) ~ radiation_therapy,
data = TCGA_BRCA_clinical_clean
)
ggsurvplot(
fit_km_radiation_therapy,
conf.int = TRUE,
data = TCGA_BRCA_clinical_clean,
pval = TRUE, #log-rank test
risk.table = TRUE, #tabella dei pazienti a rischio ad ogni intervallo di tempo
legend.title = "Radioterapia"
)
print(paste("NA pathologic_stage:", sum(is.na(TCGA_BRCA_clinical$pathologic_stage))))
## [1] "NA pathologic_stage: 22"
TCGA_BRCA_clinical_clean2 = TCGA_BRCA_clinical[!is.na(TCGA_BRCA_clinical$pathologic_stage), ]
fit_km_path_stage = survfit(
Surv(overall_survival, status) ~ pathologic_stage,
data = TCGA_BRCA_clinical_clean2
)
ggsurvplot(
fit_km_path_stage,
conf.int = TRUE,
data = TCGA_BRCA_clinical_clean2,
risk.table = TRUE,
legend.title = "Stadiazione"
)
#Valutiamo l'impatto della variabile 'radiation_therapy'
fit_cox_rad = coxph(
Surv(overall_survival, status) ~ radiation_therapy,
data = TCGA_BRCA_clinical_clean
)
summary(fit_cox_rad)
## Call:
## coxph(formula = Surv(overall_survival, status) ~ radiation_therapy,
## data = TCGA_BRCA_clinical_clean)
##
## n= 976, number of events= 110
##
## coef exp(coef) se(coef) z Pr(>|z|)
## radiation_therapyyes -0.5519 0.5758 0.1914 -2.883 0.00394 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## radiation_therapyyes 0.5758 1.737 0.3957 0.838
##
## Concordance= 0.57 (se = 0.028 )
## Likelihood ratio test= 8.32 on 1 df, p=0.004
## Wald test = 8.31 on 1 df, p=0.004
## Score (logrank) test = 8.53 on 1 df, p=0.004
Interpretazione: exp(coef) è
l’Hazard Ratio (HR).
Un HR < 1 è protettivo verso l’evento. Un HR di 0.5758 significa che
il gruppo “Radiation Therapy yes” ha un rischio di morte 0.5758 volte,
in ogni istante, rispetto al gruppo “Radiation Therapy no”.
L’RNA-seq (RNA sequencing) è una tecnologia che permette di misurare l’espressione di tutti i geni rilevabili in una cellula, in un tessuto o in un organismo.
Estrazione dell’RNA dalle cellule o dai tessuti (es. da campioni tumorali e/o sani).
Conversione in cDNA e sequenziamento (es. Illumina Hi-seq). Si ottengono milioni di “reads”.
Allineamento delle reads al genoma di riferimento.
Conteggio delle reads per gene -> si ottiene una matrice di espressione:
| Paziente 1 | Paziente 2 | Paziente 3 | … | |
|---|---|---|---|---|
| gene 1 | 456 | 38 | 35 | |
| gene 2 | 12 | 11 | 77 | |
| … |
#Query per i dati di RNA-seq sequenziati con Illumina HiSeq
TCGA_BRCA_RNA = getLinkedOmicsData(
project = "TCGA-BRCA",
dataset = "RNAseq (HiSeq, Gene level)"
)
#Manipolazione del dataset per avere il dataframe campioni x geni
TCGA_BRCA_RNA = t(TCGA_BRCA_RNA) %>%
as.data.frame() %>%
janitor::row_to_names(row_number = 1) %>%
mutate(across(everything(), ~ type.convert(as.character(.), as.is = TRUE)))
dim(TCGA_BRCA_RNA)
## [1] 1093 20155
# Selezioniamo solo i pazienti a cui sono associati dati di sopravvivenza
matching_ids = intersect(rownames(TCGA_BRCA_clinical), rownames(TCGA_BRCA_RNA))
BRCA_RNA = TCGA_BRCA_RNA[matching_ids,]
BRCA_clin = TCGA_BRCA_clinical[matching_ids,]
Genericamente prima di effettuare un’analisi multivariata si riducono i predittori con uno screening delle variabili in modo da avere \(d < p\) predittori informativi.
Selezioniamo alcuni geni che dalla letteratura sappiamo essere markers per il Breast Cancer: Breast Cancer gene 1 (BRCA1), Breast Cancer gene 2 (BRCA2), Estrogen Receptor Alpha (ESR1), Progesterone Receptor (PGR), Human Epidermal Growth Factor Receptor 2 (ERBB2)
gene_sign = c("BRCA1", "BRCA2", "ESR1", "PGR", "ERBB2")
BRCA_sign = BRCA_RNA[,gene_sign]
dim(BRCA_sign)
## [1] 1051 5
fit_cox_sign = coxph(
Surv(BRCA_clin$overall_survival, BRCA_clin$status) ~
BRCA1 + BRCA2 + ESR1 + PGR + ERBB2,
data = BRCA_sign
)
summary(fit_cox_sign)
## Call:
## coxph(formula = Surv(BRCA_clin$overall_survival, BRCA_clin$status) ~
## BRCA1 + BRCA2 + ESR1 + PGR + ERBB2, data = BRCA_sign)
##
## n= 1051, number of events= 151
##
## coef exp(coef) se(coef) z Pr(>|z|)
## BRCA1 0.120805 1.128404 0.086763 1.392 0.1638
## BRCA2 0.019034 1.019216 0.077641 0.245 0.8063
## ESR1 0.046387 1.047480 0.036998 1.254 0.2099
## PGR -0.069020 0.933308 0.033153 -2.082 0.0374 *
## ERBB2 0.006612 1.006634 0.054731 0.121 0.9038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## BRCA1 1.1284 0.8862 0.9519 1.338
## BRCA2 1.0192 0.9811 0.8753 1.187
## ESR1 1.0475 0.9547 0.9742 1.126
## PGR 0.9333 1.0715 0.8746 0.996
## ERBB2 1.0066 0.9934 0.9042 1.121
##
## Concordance= 0.592 (se = 0.028 )
## Likelihood ratio test= 8.27 on 5 df, p=0.1
## Wald test = 8.47 on 5 df, p=0.1
## Score (logrank) test = 8.51 on 5 df, p=0.1
Da questa analisi evinciamo che l’unico gene con p-value inferiore alla soglia convenzionale di \(\alpha=0.05\), è PGR. La sua espressione è statisticamente associata alla sopravvivenza globale.
In particolare, guardando a \(exp(coef)\), per ogni aumento unitario nell’espressione di PGR, il rischio di morte viene moltiplicato per 0.9333, mantenute costanti le altre variabili.
L’indice prognostico è il predittore lineare PI = (coef_gene1 * expr_gene1) + (coef_gene2 * expr_gene2)…
merged_data = cbind(BRCA_sign, BRCA_clin)
prognostic_index = predict(
fit_cox_sign,
type = "lp", #linear prediction
newdata = merged_data
)
# Stratifichiamo i pazienti in base alla mediana del PI
merged_data$risk_group = ifelse(
prognostic_index > median(prognostic_index, na.rm=TRUE),
"Alto Rischio",
"Basso Rischio"
)
fit_km_risk = survfit(Surv(overall_survival, status) ~ risk_group, data = merged_data)
ggsurvplot(fit_km_risk,conf.int = TRUE, pval = TRUE, title = "Stratificazione per gruppi di rischio (Cox-Multivariata 5 Geni)")
Perché glmnet?
Problema (p > n): Quando abbiamo più predittori (geni, \(p\)) che osservazioni (pazienti, \(n\)), la regressione di Cox standard fallisce.
Soluzione: Regressione Penalizzata.
glmnet implementa vari tipi di penalizzazioni, tra cui:
\[ \hat{\beta} = \arg\max_{\beta} \left\{ \ell(\beta) - \lambda P(\beta) \right\} \]
\[ P(\beta) = \alpha \|\beta\|_1 + (1-\alpha)\|\beta\|_2^2 \]
Da uno screening delle variabili abbiamo identificato i geni più variabili (i.e. potenzialmente più informativi).
#Selezioniamo HV_genes (Highly Variable Genes)
HV_genes = read.csv("HV_genes.csv")$x
BRCA_HV = BRCA_RNA[, HV_genes]
# glmnet richiede una matrice per x e un oggetto Surv per y
x_matrix = as.matrix(BRCA_HV)
y_surv = Surv(BRCA_clin$overall_survival, BRCA_clin$status)
#glmnet applica LASSO impostando $alpha=1$
cvfit_lasso = cv.glmnet(
x_matrix, y_surv,
family = "cox", alpha = 1,
nfolds = 5, # cv.glmnet trova il lambda ottimale con Cross-Validation (5-folds)
type.measure = "C" # usiamo come metrica di valutazione il C-index
)
plot(cvfit_lasso)
#Variabili selezionate
coefficients_lasso = coef(cvfit_lasso, s="lambda.min")
selected_features = rownames(coefficients_lasso)[coefficients_lasso[,1] != 0]
length(selected_features)
## [1] 144
print(selected_features[1:20])
## [1] "ACTL8" "ADAM20" "ADAMTS8" "AKR1E2" "ALS2CR12" "AMH"
## [7] "ASAH2" "BAGE2" "BCHE" "BCL2L10" "BCL8" "C14orf50"
## [13] "C1orf168" "C1orf194" "C2orf39" "C2orf82" "C3orf45" "C6orf142"
## [19] "C6orf81" "C7orf54"
glmnet non ha un sistema interno per il tuning di \(alpha\) nel caso dell’ EN, quindi è
possibile scrivere un ciclo per testare diversi alpha tra 0 e 1.
#alpha = seq(0.1, 1, length = 10)
#fitEN = list()
#for (i in 1:length(alpha)) {
#fitEN[[i]] = cv.glmnet(x_matrix, y_surv, family = "cox", alpha = alpha[i], nfolds = 5, type.measure = "C")
#}
# Troviamo il modello (alpha) con il C-index medio (cvm) migliore
#idx = which.max(sapply(fitEN, function(xx) {
#xx$cvm[xx$lambda == xx$lambda.min]
#}))
#cvfit_EN = fitEN[[idx]]
#best_alpha = alpha[idx]
Per semplificare impostiamo \(alpha=0.5\)
cvfit_EN = cv.glmnet(
x_matrix, y_surv,
family="cox", alpha=0.5,
nfolds=5, type.measure="C"
)
plot(cvfit_EN)
coefficients_EN = coef(cvfit_EN, s="lambda.min")
selected_features_EN = rownames(coefficients_EN)[coefficients_EN[,1] != 0]
length(selected_features_EN)
## [1] 609
print(selected_features_EN[1:30])
## [1] "ABCA4" "ABCC6P1" "ABCG4" "ABP1" "ACRV1" "ACSBG2"
## [7] "ACTA1" "ACTBL2" "ACTC1" "ACTL8" "ACTN3" "ADAM20"
## [13] "ADAMTS19" "ADAMTS8" "ADRA1A" "AKNAD1" "AKR1E2" "ALDH3A1"
## [19] "ALOXE3" "ALS2CR11" "ALS2CR12" "AMAC1L2" "AMH" "AMN"
## [25] "AMY1A" "ANKRD2" "ANO3" "ANO5" "ANXA13" "APOF"
Calcoliamo i Prognostic Index (usando lambda.min, il lambda ottimale)
lp_lasso = predict(cvfit_lasso, newx=x_matrix, s="lambda.min", type="link")
lp_en = predict(cvfit_EN, newx=x_matrix, s="lambda.min", type="link")
# Stratifichiamo i gruppi di rischio
merged_data$risk_group_lasso = ifelse(lp_lasso > median(lp_lasso), "High", "Low")
merged_data$risk_group_enet = ifelse(lp_en > median(lp_en), "High", "Low")
fit_km_lasso = survfit(Surv(overall_survival, status) ~ risk_group_lasso, data=merged_data)
ggsurvplot(fit_km_lasso, conf.int = TRUE, pval=TRUE, risk.table = TRUE, title="LASSO (alpha=1)", legend="none")
fit_km_enet = survfit(Surv(overall_survival, status) ~ risk_group_enet, data=merged_data)
ggsurvplot(fit_km_enet, conf.int = TRUE, pval=TRUE, risk.table = TRUE, title="Elastic Net (alpha=0.5)", legend="none")
Grazie per l’attenzione!
Per approfondimenti e maggiori informazioni su metodi avanzati per l’analisi di sopravvivenza siete invitati alla Poster Session: Bridging Multi-omics Data ad Survival Analysis with Omics2Surv.
Filtra i pazienti con dati mancanti per la variabile
“PAM50”.
Fitta il modello.
Visualizza la kaplan meier stratificata.
Interpreta il log-rank test.
Includi la variabile “radiation_therapy” alla cox
multivariata.
Interpreta i risultati.
Stratifica per gruppi di rischio.
Prova a lanciare il codice per il tuning di alpha.
This tutorial is part of the dissemination activities supported by the P2022BLN38 project Computational Approaches for the integration of multi-omics data – funded by European Union – Next Generation EU within the PRIN 2022 PNRR program (D.D. 1409 del 14-09-2022 Ministero dell’Università e della Ricerca) CUP B53D23027810001.